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cal  solution  of  the  shallow  water  equations  along  with  the  stability  of  these  methods. 
Most  of  the  thesis  is  concerned  with  the  background  and  formulation  of  the  shallow 
water  equations.  The  derivation  of  the  basic  equations  will  be  given,  in  the  primitive 
variable  and  vorticity-divergence  formulation.  Also  the  shallow  water  equations  will 
be  written  in  spherical  coordinates.  Two  main  types  of  methods  used  in  approximat¬ 
ing  differential  equations  of  this  nature  will  be  discussed.  The  two  schemes  are  finite 
difference  method  (FDM)  and  the  finite  element  method  (FEM).  After  presenting 
the  shallow  water  equations  in  several  formulations,  some  examples  will  be  presented. 
The  use  of  the  Fourier  transform  to  find  the  solution  of  a  semidiscrete  analog  of  the 
shallow  water  equations  is  also  demonstrated. 
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I. 


INTRODUCTION 


The  majority  of  this  work  is  a  compilation  of  past  mathematical  papers  con¬ 
cerning  the  numerical  solution  of  the  shallow  water  equations,  and  I  can  in  no  way 
take  credit  for  them.  Hopefully  I  have  put  this  material  in  a  readable  and  under¬ 
standable  text  that  will  provide  a  good  basis  for  others  to  use  for  reference  in  this 
area  of  mathematics.  The  next  several  paragraphs  will  highlight  the  key  points  in 
slightly  more  detail  than  the  Table  of  Contents. 

Chapter  II  shall  review  and  discuss  the  mathematical  formulation  of  flow 
in  shallow  regions.  Two  formulations  will  be  given,  namely  the  primitive  and  the 
vortirity-divergence  formulations  of  the  shallow  water  equations.  It  will  also  look  at 
the  formulation  of  the  equations  in  spherical  coordinates. 

Cljapter  III  will  look  at  the  linearization  of  the  shallow  water  equations.  In 
Cartesian  coordinates  we  discuss  the  one  dimensional  and  two  dimensional  cases. 
Tlie  linearized  vorticity-divergence  form  of  the  shallow  water  equations  is  derived  in 
this  chapter  for  the  two  dimensional  case.  The  spherical  case  is  also  discussed  here. 
Note  that  in  the  Cartesian  case,  the  linearized  equations  have  constant  coefficients  in 
contrast  to  the  spherical  coordinates.  This  means  that  linear  stability  analysis  will 
be  more  difficult  in  the  latter. 

Chapter  IV  shall  develop  two  of  the  different  types  of  approximations  com¬ 
monly  used.  The  first  is  the  finite  difference,  and  the  second  is  the  finite  element.  The 
finite  element  may  also  be  referenced  to  as  the  Galerkin  method  and  will  be  referred 
to  l>y  both  names  in  this  thesis.  Other  families  of  methods  such  as  spectral  or  finite 
volume  approximation  are  beyond  the  scope  of  this  thesis. 

In  C'liapter  V,  the  stability  analysis  of  the  shallow  water  equations  will  be 
examined.  For  this  analysis  we  use  the  linearized  equations  obtained  in  C’liapter  II. 
Fourier  transforms  in  a  one  dimensional  case  will  be  covered.  Fourier  transforms 
will  also  be  used  to  look  at  the  two  dimensional  case.  The  spherical  case  requires 
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special  consideration  as  in  Longuet-Higgins(see  [Ref.  1])  or  Neta  (see  [Ref.  2]).  A 
summary  of  the  results  for  the  various  techniques  are  provided  in  a  table  for  ease  of 
later  examination. 

Also  covered  in  this  report  will  be  an  extensive  bibliography  and  reference 
section,  so  that  future  students  studying  in  this  area  will  be  able  to  pick  up  more 
easily  where  this  thesis  leaves  oif. 
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II. 


SHALLOW  WATER  MODEL 


A.  MODEL  BACKGROUND 

For  this  model  consider  a  sheet  of  fluid  with  constant  and  uniform  density. 
(See  [Ref.  3].)  The  height  of  the  surface  of  the  fluid  above  the  reference  level  z  =  0  is 
h{x,  y,  t).  We  model  the  body  force  arising  from  the  potential  (p  =  gh  with  atmosphere 
or  ocean  in  mind.  ^  is  a  vector  directed  perpendicular  to  the  z  =  0  surface,  or  (j)  can 
be  said  to  be  antiparallel  to  the  vertical  axis  (i.e.,  <!>  is  in  the  direction  opposite  to 
the  vertical  axis).  The  rotation  axis  of  the  fluid  is  the  z-axis  in  this  model.  In  this 
ease  the  Coriolis  parameter  /  is  217  since  J7  =  fcfl.  The  rigid  bottom  is  defined  by 
the  surface  z  =  /is(x,y).  The  velocity  has  components  u,  v,  and  w  in  the  x,  y,  and 
z  directions  respectively.  Though  the  pressure  of  the  fluid  surface  can  be  arbitrarily 
imposed,  for  this  model  it  will  be  assumed  to  be  constant.  Lastly,  the  fluid  is  assumed 
inviscid,  in  other  words,  only  the  motions  for  which  viscosity  is  not  important  are 
considered. 

In  this  model,  because  the  depth  of  the  fluid,  h  —  hs,  varies  over  time  or  space, 
let  H  be  the  average  depth  of  the  fluid.  H  characterizes  the  vertical  scale  of  motion 
also.  Let  L  be  the  characteristic  horizontal  scale  for  the  motion.  Then  a  fundamental 
condition  which  will  characterize  shallow  water  theory  will  be 


which  is  also  called  the  hydrostatic  approximation  with  long  wavelengths.  (See  [Ref. 
4].)  The  shallow  water  model  contains  several  of  the  important  dynamical  features 
of  the  atmosphere  and  ocean  while  being  simple  enough  to  be  easily  understood. 
The  major  physical  difference  with  this  model  and  reality  is  the  absence  of  density 
stratification  that  is  present  in  the  real  fluids  such  as  Earth’s  atmosphere  or  oceans. 
The  hydrostatic  approximation  also  allows  p  to  vary  with  z,  but  we  will  consider  p 
constant  for  this  model. 
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Recalling  the  equation  for  motions  of  fluids  from  Haltiner  and  Williams  (see 
[Ref.  5],  we  have 

^  =  -aVj9-2fi  X  +  (2.1) 

o  is  the  specific  gravity  (i.e.  ^),  g  is  the  sum  of  gravitational  and  centrifugal  forces, 

^  — * 

and  F  is  the  force  due  to  friction.  For  our  model  we  will  be  neglecting  g  and  F 

because  with  our  assumptions  they  are  much  smaller  in  magnitude  than  the  Coriolis 

forces.  Recall  that 

dV  dV  I  ^  ^  -  - 

^  =  ^  +  5V(V.r)  +  (VxV)xV.  (2.2) 

^  w  w  — * 

where  V''  =  ui  +  vj  +  wk.  Using  these  two  basic  equations,  we  will  construct  our 
model. 

B.  MODEL  EQUATIONS 

There  are  several  formulations  in  the  literature.  We  will  cover 

•  primitive  variable  formulation 

•  vorticity  divergence  formulation. 

The  formulations  using  stream  function  and  velocity  potential  will  not  be  discussed 
here,  (see,  for  example,  [Ref.  5]). 

1.  Primitive  Form 

Now  follow  the  consequences  of  the  model  in  the  realm  of  the  dynamical  equa¬ 
tions  of  motion.  Recalling  the  Navier-Stokes  equations  which  describe  tin’  conser¬ 
vation  of  mass  and  momentum  (see  [Ref.  6]),  we  see  that  the  dynamics  and  the 
thermodynamics  decouple  due  to  the  specification  of  incompressibility  and  constant 
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density.  (See  [Ref.  3].)  This  reduces  the  equation  of  mass  conservation  to  the  condi¬ 
tion  of  incompressibility: 

fi'it  fi't) 

(2.3) 


du  dv  ^ 

dx^  dy^  dz 


U 


Now  the  first  two  terms  of  this  equation  are  of  order  where  U  can  be  considered 
the  characteristic  scale  for  the  horizontal  velocity.  It  then  follows  that  the  scale  for 
the  vertical  velocity  (W^)  is  smaller  than  or  equal  to  the  order  SU.  This  represents  an 

upper  bound  for  the  vertical  velocity,  and  it  can  be  smaller  than  order  5U  if  there  is 

11-1  du  dv 

cancellation  between  and  7^. 

dx  dy 

Now  an  estimate  of  the  momentum  equations  in  component  form  are: 
du  du  dudu  1  dp 

^  +  +  +  =  —pai 

U  f/2  f/2  UW  P 

T  L  L  H  ^  pL 


dv  dv  dv  dv  .  1  dp 

+  +  =  —pd-y 

U  {/2  U2  p 

T  L  L  H  pL 


(2.4) 


dw 

dt 

W 


dw  dw 
+  U—  -f  v—  +  w 
dx  dy 

UW  UW 


dw 

dz 

W^ 


1  dp 
p  dz 

P 


T  L  L  H 


pH 


f  will  be  defined  differently  depending  on  the  form  of  the  shallow  water  equations. 


Ill  (2.4)  each  term  has  the  order  of  magnitude  written  immediately  below  it  in  terms 
of  the  characteristic  scales  where  T  is  the  characteristic  scale  for  time  and  P  is  the 


characteristic  scale  for  the  pressure  field.  Equation  (2.4)  follows  from  (2.1)  and  (2.2) 


ill  the  following  manner.  Remember  that  we  are  neglecting  the  gravity  term  and  the 
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force  due  to  friction.  Now  when  we  expand  out  (2.1),  we  get 

dV  \dp-f  \dp-t 

dt  pdx^  p  dy^  p  dz 


+20(j  cos  c/7  +  sin  (/?)  x  {ui  +  vj  +  wk)^ 


which  equals 


dV  1  dp^  1  dp-t  1  dp^ 

dt  pdx^  pdy^  pdz 


(2.5) 


+{2w0,  cos  (f  —  fv)i  +  fuj  —  2uQ,  cos  (pk. 

Since  we  are  modeling  a  rectangular  system  first,  the  p  terms  drop  out  and  (2.5) 
becomes 


dV  1  dp^  Idp-t  Idpp  r  ^  ,  r  lo  a\ 

pox  pOy  poz 


dt 


When  we  look  at  the  spherical  form  of  the  shallow  water  equations,  the  <p  terms  are 
used.  Note  that 

-  -r,  ( dv\  w  / du  dw\  t  f  dv  du\  -* 


dy  dz 
Now,  expanding  out  (2.2)  we  get 


dt 


dV 


du  dv  dw\  -f 


,dx  dy^ 

du  dv  _  dw^ 


—  -  —+  M— +  h  + 


dt 


dx  dx 


dx 


dy  ■  dy  dy, 


+ 


^  du  dv 
u-^  + 

i  uz  Oz 


k  + 


[dz 


dw  (  dv 


i 


Collecting  terms,  this  becomes 


dt 


du  du  du  du\ 
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or  more  simply 


dt 


Equating  the  right  hand  side  of  (2.6)  and  this  equation  and  then  separating  into  the 
three  components  will  give  us  (2.4). 

Note  that  the  total  pressure  (p)  is 


p{x,  y,  z,  t)  =  -pgz  +  p{x,  y,  z,  t). 


The  horizontal  pressure  gradient  is  independent  of  z  if  one  uses  the  hydrostatic  ap¬ 
proximation 

Y^  =  -pg  +  o(i^)  (2.7) 

where  (J  is  —  and  J  1  .  This  approximation  follows  from  the  scale  analysis  of 
L/ 

the  momentum  equations  and  the  incompressibility  condition.  Integrating  (2.7)  and 
using  the  boundary  condition 

p(x,  y,  h)  =  po 


yields 


p  =  pg{h  -z)+  Po. 


A  way  to  look  at  this  is  to  think  of  the  pressure  which  is  in  excess  of  po  at  any  point 
simply  as  the  weight  of  the  unit  column  of  fluid  above  that  point  at  that  particular 
instant  in  time.  Remember  that  the  horizontal  pressure  gradient  is  independent  of  z, 
and  therefore  the  horizontal  accelerations  are  independent  of  z. 

The  Rossby-wave  number  is  the  ratio  of  inertia  to  the  Coriolis  terms  or  {^y). 
Recalling  the  Taylor-Proudman  theorem  (see,  [Ref.  3]  p.  43)  which  simply  put  states 
that  if  the  Rossby-wave  number  is  small,  friction  can  be  ignored,  and  baroclinic  vector 
is  identically  zero  (i.e.,  there  is  no  additional  pressure  change  with  a  change  in  height 
as  in  this  model),  then  it  follows  that  the  horizontal  accelerations  are  independent  of 


du  dx) 
dz  dz 
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In  this  model,  the  Taylor- Proudman  theorem  applied  to  a  homogeneous  fluid  requires 
the  velocities  to  be  independent  of  2r.  This  allows  the  horizontal  momentum  equations 


to  become 


du  dudu  1  dp 


(2.8) 


dv  dvdv  1  dp 

dt  dx  dy  p  dy 

Let  us  look  at  the  scale  analysis  of  the  vertical  momentum  equation  of  (2.4). 
Consider  w  to  be  of  O  =  0(6).  Scale  u  and  v  by  U,  and  w  by  .  By  (2.3)  we 


know  that 


where 


du  dv  « 

dx'  dy'  dz' 

f  ^ 

u  =  Uu  ^  ~  L 


V  =  Uv  y' 


y. 

L 


w 


SUw 


z'  = 


H' 


The  vertical  momentum  equation  can  be  written  as 

SU 


dw  U  dw  U . dw  ~ dw 


\  \  dp 
H  p  dz' 


Multiply  through  by  H  and  this  becomes 


W 


,,dw  rrr-^^  r  r  c  ~  CTT-^'^ 


1  dp 
p  dz' 


Removing  all  the  terms  of  0(6^),  note  from  (2.7)  that  p  is  of  O(S^)  also,  we  have 


Therefore  the  vertical  momentum  equation  is  an  identity. 


Now  we  will  rewrite  the  incompressibilty  condition  (2.3), 

dw  _  /  du 

dz  \5x  dy J 
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Since  the  horizontal  accelerations  are  independent  of  z,  we  can  integrate  the  above 
relationship  to  get 


/  -/  \ 

w{x,y,z,t)  =  -z  1^  +  ^j  -\-w{x,y,t). 

The  condition  of  no  normal  flow  at  the  bottom  requires 

/  L  j\  dhjg  dhg 
w{x,  y,  hs,  t)  =  u--  +  v--. 

ox  oy 

It  naturally  follows  that 


(2.9) 


/  S\  rr  \  ( ^'^  dhs 

v,(x,y,z,t)  =  (fc*  -  z)  +  -j  +  „—  +  . 

dz 

The  vertical  velocity  u;(x,  y,  z,t)  =  —  at  the  upper  surface  z  =  h  represents  the  rate 

at 

dh 

at  which  the  free  surface  is  rising.  Thus  w(x,y,h,t)  =  —  and  (2.9)  become 

dt 

,  ,  dh  dh  dh 


Therefore 


W  ^  Wfc  -  Ab)1  +  I;  [v(h  -  Ab)]  =  0. 

This  equation,  combined  with  the  horizontal  momentum  equations  (2.9)  form  the 
shallow  water  equations. 


du  du  du  .  dh 

a  +  “a;  +  ”5?  - 


~dt 


dv  dv 

+  u—  +  V—  +  fu 

dx  dy 


(2.10) 


()  ii  ^  d 

u,  and  h  art*  called  the  primitive  variables  and  (2.10)  is  the  primitive  form  of 
the  shallow  water  equations. 
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2.  Vorticity-divergence  Form 

For  vorticity-divergence,  we  start  with  the  definitions.  The  relative  vorticity, 
denoted  by  is  defined  by 

C  ~  ^3/5 

and  the  absolute  vorticity,  denoted  by  is  given  by 

Q  =  C  +  f.  (2.12) 

The  divergence,  denoted  by  D,  is  given  by , 

D  =  Ua,  +  Vy,  (2.13) 


and  the  kinetic  energy  (per  unit  meiss),  denoted  by  is  given  by 

^  -k 


(2.14) 


To  get  the  vorticity-divergence  formulation,  we  start  by  differentiating  the  first 
equation  of  (2.10)  with  respect  to  y  and  subtract  it  from  the  second  equation  of  (2.10) 
differentiated  with  respect  to  x  gives 


"b  UxVx  “b  y>Vxx  ”b  V^Vy  -j-  VVxy  -b  f'^x  “b  Q^xy 

Uyf  XlyXLx  DyUy  ~~  VUyy  “j”  f  Vy  Q^Xy  ~  0. 

Rearranging  and  cancelling  terms  will  yield 

(Vl  -  Uy)^  +  Ux{Vx  -  Uy)  -I-  U{Vx  -  Uy)^ 

+  Vy{Vx  -  Uy)  +  V{Vx  -  Uy)^  -j-  f  {u x  -^T  Vy)  =  0. 

Rewriting  in  terms  of  (2.11)  and  (2.13)  gives 


(2.15) 


Ct  +  +  ^Cx  +  ^yC  +  +  f  D  =  0, 

This  can  be  rewritten  in  the  form 


Ct  +  {uOx  +  i^Oy+  fD  =0 
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or 


This  is  called  the  vorticity  equation.  To  get  the  divergence  equation,  we  difFfer- 
entiate  the  first  equation  of  (2.10)  with  respect  to  x  and  add  to  it  the  second  equation 
of  (2.10)  differentiated  with  respect  to  y  gives 


Uxt  +  uj  +  UUa:x  +  WarUj,  +  VUxy  -  fv^  +  gKx 

d”  Uy'^x  ”f"  UVxy  ”1“  Vy  T  d”  f'^y  d"  ffhyy  —  0. 
Again,  rearranging  terms  will  yield 

{Ux  d-  Vy)^  +  u{ux  d-  Vy)^  +  v{Ux  d-  Vy)y-\- 

Uj,  ”1“  Vy  ”1“  ^V/yVx  fi^^x  ^y)  d”  gi^hxx  d"  hj^y)  —  0. 

Using  (2.11)  and  (2.13)  this  can  be  further  reduced  to 


(2.16) 


Dt  d-  uDx  d-  vDy  +  Ux^  d-  Vy^  +  2uyVx  —  /C  d-  =  0. 


The  vorticity-divergence  form  of  the  shallow  water  equations  in  terms  of  rela¬ 
tive  vorticity  and  divergence  are 


0  d- («C).  +  «),  d- //)  =  0 

Dt  d-  uDx  d-  vDy  +  Ux^  d-  Vy"^  2uyVx  —  /C  d-  5'V^h  =  0  (2-17) 

lit  +  [u{h  —  Hb)]^  +  [v{h  -  hB)]y  =  0. 

To  simplify  a  bit  further,  we  can  use  (2.11)  through  (2.14)  on  (2.15)  and  (2.13). 
First,  rearrange  and  cancel  terms  of  (2.15)  in  the  following  manner. 

I'xt  Hyt  d“  llxVx  UyXlx  -|-  JVx  d~  HHxx  ll'Hyx~^ 

VyVx  “  l^ylly  d“  d“  l^l^xy  HHyy  “ 

Rearranging  further  gives 

{Vx  -  Uy)^  -f  Ux{Vx  -Uy  +  f)  +  u[Vxx  -  Wyx)d- 
lly{Vx  Wy  +  f)  -|-  v{Vxy  llyy)  —  0. 
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\ 

This  leads  to 

Ct  “I”  “h  “I”  T  '^Qy  — 

or 

C.  +  {uQl  +  (vQ),  =  0.  (2.18) 

This  can  also  be  written  as 

Qt  +  "t"  ^'^Q)y  ~  ^ 

since  ft  =  0. 

Now  we  simplify  (2.16)  in  a  similar  fashion.  First,  rewrite  the  equation  in  the 
following  manner. 

"^xt  ”1"  '^yt  y^xx  ”f”  Q^yy  "f  '^^xx  ”1”  Ux'^x  "h  VVxx  ~1~  '^x‘^x~\~ 

UXiiyy  ”i"  UyUy  ~j~  Wyy  x'^ X  ”1”  ‘^x'^y  f X  '^'^XX'^ 

UUyj;  ”1"  XLy'Ox  ‘^y'^y  *1”  f  ^y  "h  '^'^xy  ~~  '^'^yy  —  ^ 

Notice  that  we  added  in  and  subtracted  some  terms  that  were  alike.  Now  reduce  this 
equation  to 

Uxt  +  Vyt  +  ghxx  +  ghyy  +  [uux  +  vVx)^-\- 
{uUy  +  VVy)y  —  VxVx  +  VxUy  “  /Vj;  “  VVxX^ 

VUyx  +  UyVx  —  UyUy  +  f Uy  +  UV xy  “  UUyy  =  0. 

Reorganize  the  terms  to  further  reduce  to 

^xt  “f  '^yt  “I"  gf^xx  d"  g^yy  d"  I  n  I  d" 

V  /  XX 

-  VxiVx  -Uy  +  f)-  v{Vx  -Uy  +  /)^d- 

yy 

Wy(jV  -Uy  +  f)+  u{Vx  -Uy  +  f)^  =  0. 

Now  usr  ihr  definitions  of  Q  and  1\  to  reduce  this  equation  to 

(^x  d”  ^y^i  d~  di^^xx  d“  ^^3/j/)  d"  ^^xx  “T  ^^yy 
"^Qx  d“  '^yQ  d“  —  0. 
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Thus  our  equation  becomes 


Dt  +  —  {vQ)x  +  {uQ)y  =  0.  (2.19) 

Combining  (2.18),  (2.19),  and  the  last  equation  of  (2.10)  we  have  the  shallow 
water  equations  in  the  vorticity  divergence  form. 

C<  +  («Q)a;  +  {‘>jQ)y 

Dt  +  gV^h  +  V^K  -  (vQ)^  +  {uQ)y 
ht  +  [u{h  -  Hb)]^  +  [v{h  -  hB)]y 


=  0 

=  0  (2.20) 
=  0 


C.  SPHERICAL  COORDINATES 
1.  Background 

Where  Cartesian  coordinates  locate  a  point  using  3  vectors,  spherical  coor¬ 
dinates  locate  the  same  point  using  two  angles  and  a  distance.  Typically,  spherical 
coordinates  represent  a  point  in  space  with  an  ordered  triplet  such  as  (a,  <^,  A).  These 
variables  are  related  to  (x,y,z)  by 

X  =  a  sin  ^  cos  A 

y  =  a  s'm  (f>  sm  X  (2.21) 

z  =  a  cos  (j). 

Let  us  consider  a  channel  that  goes  around  the  whole  Earth.  In  this  model, 
let  A  be  the  longitude,  let  cj)  be  the  latitude,  let  r  be  the  radial  distance,  let  a  be  the 
average  radius  of  the  Earth,  and  let  z  be  the  average  sea  level.  Now,  recall  equations 
(2.5)  and  (2.2). 

dV  _  \  dp-f  \dp-t 

dt  pdx  pdy^  pdz 

-f  (2tnf2  cos  d)  —  fv)i  fuj  —  2uQ,  cos  (pk. 
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and 


^  =  f +  iv(v.v)  +  (vxy)xp. 

Using  results  drawn  from  Haltiner  and  Williams  (see  [Ref.  5]),  and  neglecting  the 
force  due  to  friction,  F,  gives  us 


f  du  du  du  du\  ^ 


( dv  dv  dv  dv\  ->  ( dw  dw  dw  dw\ 


u 


+  {i|?+(2(l+  ^ 

ypax  \  a  cos  (pj 


(v  sincf)  —  w  cos  <l>)\  i 


(2.22) 


+  2n  +  -^ 

(poy  \  a  cos  (pj 


vw  -t 

usm<p - >  j 

a 


+ 


+  5^  d"  1  + 


u  \  ,  ^  r 

- -  u  cos  (p-\ - >  k 

a  cos  (j)j  a  \ 


Separating  along  the  three  axes  and  we  have 

.  du  du  du  du  1  dp  u  '' 

0  =  —  +  u— — I-  V— — h  w— — I — — - j  2f2  H - 7 

dt  dx  dy  dz  p  dx  \  a  cos  <p^ 

^  dv  dv  dv  dv  1  dp  u  \ 

0  =  —  +  u— — |-  u-7 — h  w— — I — -7 — h  I  2f2  H - 7  1 

dt  dx  dy  dz  pdy  \  a  cos  <p  J 


[v  sin  (f)  —  w  cos  <^) 


vw 


usmo  — 


(2.23) 


„  dtr  dw  dw  dw  1  dp  /  ^  u  \  , 

0  =  —  +  u—  +  V—  +  w-^  +  -\-g+  2U  + - 7  I  ucos</)+  — 

dt  dx  dy  dz  pdz  ^  a  cps  q>J  a 

.  T  ,  .  .  I  1  r  J  1  1  1  1  .  .  r  dh 

Note  that  p  is  independent  of  z.  and  and  can  be  written  m  terms  of  p— 

pdx  pdy  dx 

^  •  .  A. 

and  g—  respectively.  Also  note  that 
dy 

du  du  dX  du  d(j) 

dx  dX  dx  ^  d4>dx 


du 


1 


+ 


du 


1 


5A  — asind>sin  A  d<f>  a  cos  (f)  cos  X 
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and  similarly  for  the  other  partials.  Substituting  these  and  the  Coriolis  parameter, 
/  =  2fisin0,  into  (2.23)  gives  us  the  shallow  water  equations  in  spherical  coordinates. 
(See  [Ref.  7].) 


du  1 
dt  a  cos  0 


+  V  cos  6 


du 

'm 


/  +  —  tan  6 
a 


9  dh 
a  cos  6  dX 


dv  ^  1 

dt  a  cos  6 


dv  dv 


+ 


u 


/  H —  tan  0 
a 


9  dh 
a  dd 


(2.24) 


dt  a  cos  6 


+  e) 


=  0. 


Numerical  approximations  of  the  shallow  water  equations  in  spherical  coordinates  are 
given  for  example  by  Turkel  and  Zwas  (see  [Ref.  8])  and  Arakawa  (see  [Ref.  9])  et 
al.  Analysis  of  the  schemes  were  given  by  Neta  and  Navon  (see  [Ref.  10])  et  al. 
Navon  and  de  Villiers  et  al.  have  applied  the  Turkel-Zwas  scheme  to  a  hemispheric 
barotropic  model  (see  [Ref.  11]). 
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Ill 


LINEARIZATION 


In  this  chapter,  we  will  see  how  to  linearize  the  shallow  water  equations  in 
Cartesian  and  spherical  coordinates.  The  Cartesian  coordinate  case  starts  with  either 
the  primitive  or  the  vorticity-divergence  formulation. 

A.  TWO  DIMENSIONS 
1.  Basic  Formula 

First  separate  each  of  the  three  variables  into  the  average  which  we  can  consider 
constant  plus  the  perturbation  in  the  variable.  This  should  look  like: 

u  =  U  +  u' 

V  =  V  +  v'  (3-1) 

h  =  H  +  h' 

where  the  capital  letters  are  used  for  the  average  or  mean  values  and  the  primed 
variables  are  the  perturbations.  Now  substituting  this  into  (2.7)  and  omitting  all  the 
second  order  terms  gives 


du' 

dt 


du'  du' 


fV-fv'  +  g 


d{H  +  h') 
dx 


0 


i)v'  ,.dv' 


,dv' 


d(H  +  h') 
dy 


=  0 


(3.2) 

^(>r  simplicity  and  without  loss  of  generality  we  sometimes  consider  a  zero 
mean  flow,  i.e.  I  =  V""  =  0.  This  and  dropping  the  primes  will  give  the  linearized 
form  of  the  shallow  water  equations  which  follows. 
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du  dh 


=  0 


dv 


z-  dh  ^ 
+  fu  +  g-^  =  0 
oy 


dh 

j-^+(H-hB) 


/ du  _ 

\5x  dy) 


0. 


If  topography  is  not  considered,  then  the  linearized  shallow  water  equations 
can  be  written 

du  dh 


dv  dh 


This  can  also  be  written  as 


dh  jj  (d'^ 

dt  ^  ^  dyj 


=  0. 


ut  —  fv  +  ghx  =  0 

4-  /u  +  ghy  —  0  (3.3) 

ht  +  H{ux  +  Vy)  =  0. 

Another  source  for  the  development  of  the  linearized  shallow  water  equations  is  Gill. 
(See  [Ref.  4]). 

2.  Vorticity  Divergence  Formula 

The  two  dimensional  linearized  vorticity  divergence  form  of  the  shallow  water 
equations  is  obtained  by  using  (2.8)  and  (2.9)  on  (3.3).  It  should  be  noted  that 
the  linearized  vorticity  divergence  formula  can  be  found  linearizing  (2.16)  or  (2.19). 
Taking  the  partial  differential  of  the  first  equation  of  (3.3)  with  respect  to  y  and 
subtracting  it  from  the  partial  differential  of  the  second  equation  of  (3.3)  with  respect 
to  X  yields 

Vtx  +  fUx  +  ghxy  -  Uty  +  fVy  -  ghxy  =  0. 
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With  a  little  algebra  this  equation  can  be  rewritten  as 

(Ua:  -  Uy)^  -I-  f{Ua:  +  V^c)  =  0. 

which  when  using  the  definitions  of  (  and  D  reduces  to 

Ct  +  fD  =  0.  (3.4) 

This  is  the  first  equation  of  the  linearized  vorticity  divergence  form.  By  taking  the 
partial  of  the  first  equation  of  (3.3)  with  respect  to  x  and  adding  it  to  the  partial 
differential  of  the  second  equation  of  (3.3)  with  respect  to  y  we  get 

utx  -  fvx  +  ghxx  +  vty  +  fuy  +  ghyy  =  0. 

With  a  little  algebraic  simplification  our  equation  becomes 

(ux  T  fiy^x  ^y)  T  gij^xx  T  ^yy)  “ 

Again,  using  the  definitions  of  C,  and  D,  this  reduces  to 

A  -  /C  +  g'^^h  =  0.  (3.5) 

(3.4)  and  (3.5)  along  with  the  last  equation  of  (3.3)  form  the  two  dimensional  lin¬ 
earized  vorticity  divergence  form  of  the  shallow  water  equations. 

0  +  /T>  =  0 

A-/C  +  ^VA  =  0  (3.6) 

ht  +  H(ux  +  Vy)  =  0 
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B.  LINEARIZED  SHALLOW  WATER  EQUATIONS  IN 
ONE  DIMENSION 


with 


In  the  one  dimensional  case,  all  dependence  on  y  is  eliminated  and  we  end  up 


du 


dv 

'm 


+  fu  =  0 


or  more  simply 


=  0, 


ut-  fv+  gK  =  0 
vt  +  fu  =0 

hf  —  0. 


(3.7) 


C.  SPHERICAL 

To  linearize  the  spherical  version  of  the  shallow  water  equations,  simply  take 

(3.1)  and  substitute  it  into  (2.26).  The  resulting  equation  follows  quite  easily  once 

all  second  and  higher  order  terms  are  omitted. 

du  gdh 

H - —  0 

ot  a  cos  B  aX 


dv 

m 


r  9  dh  ^ 


m  H 

dt  ^  a  cos  0  \dX 


(3 


.Notice  that  the  linearized  equations  have  nonconstant  coefficients.  This  complicates 
the  linear  stabilty  analysis.  See  Longuet-Higgins  [Ref.  1]  and  Neta  [Kef.  12]  et  al. 
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IV.  APPROXIMATIONS 


A.  FINITE  DIFFERENCES 

1.  Introduction 

One  of  the  first  steps  in  using  finite  difference  methods  is  to  replace  the  con¬ 
tinuous  problem  domain  by  a  difference  mesh  or  a  grid.  Let  f{x)  he  a  function  of 
the  single  independent  variable  a:  for  a  <  a;  <  6.  The  interval  [a,  6]  is  discretized  by 
considering  the  nodes  a  =  xo  <  xi  <  •  ■  •  <  xn  <  a^jv+i  =  and  we  denote  f{xi)  by 
The  mesh  size  is  Xj+i  —  Xj  and  we  shall  assume  for  simplicity  that  the  mesh  size 
is  a  constant 

,  b  —  a 

'■“/vTT 

and 

Xi  =  a  ik  f  =  0, 1 ,  •  •  • ,  -|- 1 

In  the  two  dimensional  case,  the  function  f{x^y)  may  be  specified  at  nodal 
point  {xi,yj)  by  fij.  The  spacing  in  the  x  direction  is  hx  and  in  the  y  direction  is  hy. 

Taylor  series  expansions  of  functions  in  several  variables  play  an  important 
role  in  formulation  and  classification  of  finite  difference  methods.  The  Taylor  series 
expansion  for  /,+!  about  the  point  x,-  is  given  by 

fi+i  =  /»■  +  - 

The  Taylor  series  expansion  for  about  the  point  {xi,yj)  is  given  by 

/l2  /j2 

/ii+l  j+I  ~  fij  d"  i^^xfx  d"  ^yfy)ij  d"  A®  d"  hxhyfxy  H  ~ fyy)ij  d"  '  ‘  ‘ 

2.  Finite  Differences 

Using  text  drawn  from  Neta  (see  [Ref.  13]),  we  will  cover  certain  aspects  of 
finite  differencing.  An  infinite  number  of  difference  representations  can  be  found  for 
the  partial  derivatives  of  /(x,t/).  Let  us  use  the  following  operators: 

Ax/tj  =  /i+i/  —  fij  1*^  forward  difference  operator 
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Vxfij  =  fij  —  fi-ij  1®*  backward  difference  operator 
Sxfij  =  fi+ij  —  fi-ij  centered  difference 
^xfij  —  fi+l/2j  ~ 

=  (/i+i/2i  +  /i-i/2i)/2  averaging  operator 

Note  that 

In  a  similar  fashion  we  can  define  the  corresponding  operators  in  y. 

In  the  following  table  we  collected  some  of  the  common  approximations  for 
the  first  derivative. 


Table  I.  First  Derivative  Approximations 


Finite  Difference  Order 

(See  Chapter  IV  Section  A  3) 

0{hx) 

fix 

0(h.) 

0(hl) 

+  4/,+„  -  /i+y)  =  0(hl) 

5^(3/.,  -  4/,-„  +  =  i(V.  +  iv=)/„  0.(kl) 

0(hl) 


The  compact  fourth, order  three  point  scheme  deserves  some  explanation.  Let 
fx  be  ?’,  then  the  method  is  to  be  interpreted  as 

(1  + 
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or 

g  H"  4:Vij  +  * 

df 

This  is  an  implicit  formula  for  the  derivative  —  at  {xi^  yj).  The  Vij  can  be  computed 

ox 

from  the  fij  by  solving  a  tridiagonal  system  of  algebraic  equations. 

The  most  common  second  derivative  approximations  are 

1 


fxx\ij  ^2  ”1“  fi+2j)  4“  0[hx) 


fxx\ij  L 2  “f  /t-2j)  "h  0(^hx) 

^x 

fxx\ij  = 


fxx\ij  — 


1  Slfn 


hi  1  +  i-^si 


+  0(hl) 
+  0(hi). 


Remarks: 

1.  The  order  of  a  scheme  is  given  for  a  uniform  mesh. 

2.  Tables  for  difference  approximations  using  more  than  three  points  and 
approximations  of  mixed  derivatives  are  given  in  Anderson,  Tannehill  and  Fletcher 
(see  [Ref.  14]. 


3.  Difference  Representations  of  PDEs 


a.  Truncation  Error 

The  difference  approximations  for  the  derivatives  can  be  expanded  in 
Taylor  series.  The  truncation  error  (T.E.)  is  the  difference  between  the  partial  deriva¬ 
tive  and  its  finite  difference  representation.  For  example 


/x 


-  =  /. 

ij  hx 


u 


/t-flj  fij  _  _  f 


hx 


We  use  0{lix)  which  means  that  the  truncation  error  satisfies  |T.  E.\  <  K\hx\  for 
hr  0,  sufficiently  small,  where  I\  is  a  positive  real  constant.  Note  that  0{hx)  does 
not  tell  us  the  exact  size  of  the  truncation  error.  If  another  approximation  has  a 
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truncation  error  of  0{h?^),  we  might  expect  that  this  would  be  smaller  only  if  the 
mesh  is  sufficiently  fine. 

We  define  the  order  of  a  method  as  the  lowest  power  of  the  mesh  size  in 
the  truncation  error.  Thus  Table  1  gives  first  through  fourth  order  approximations. 

The  truncation  error  for  a  finite  difference  approximation  of  a  given 
PDE  is  defined  as  the  difference  between  the  two.  For  example,  if  we  approximate 
the  advection  equation 


dF  dF  ^ 

^  +  c^-0,  c>0 


by  centered  differences 


Fij+i  Fij-\ 


2At 


then  the  truncation  error  is 


^  ^  (dF  dF\ 


2Ax 

_ Fjj+i  Fjj—i 

2At 


Fi^ij  -  Fi- 


—  c- 


2Ax 


\  d^F  \  d^F  . 

=  — At^-^ —  c-Aa:^-— —  —  higher  powers  of  At  and  Ax 
6  ot'^  6  ox-^ 

.  We  will  write 

T.E.  =  0{AF,Ax^). 
b.  Consistency 

A  difference  equation  is  said  to  be  consistent  or  compatible  with  the 
partial  differential  equation  when  it  approaches  the  latter  as  the  mesh  sizes  approaches 
zero.  This  is  equivalent  to 

T.E.  0  as  mesh  sizes  -f  0 . 


c.  Stability 

.\  numerical  scheme  is  called  stable  if  errors  from  any  source  (e.  g.  trun¬ 
cation.  rtnind-off.  errors  in  mecisurements)  are  not  permitted  to  grow  cis  the  calculation 
proceeds.  Hirhtnieyer  and  Morton  give  a  less  stringent  definition  of  stability  (see  [Ref. 
15]).  A  scheme  is  stable  if  its  solution  remains  a  uniformly  bounded  function  of  the 
initial  .state  for  all  sufficiently  small  At. 


d.  Convergence 

A  scheme  is  called  convergent  if  the  solution  to  the  finite  difference 
equation  approaches  the  exact  solution  to  the  PDE  with  the  same  initial  and  boundary 
conditions  as  the  mesh  sizes  apporach  zero.  Lax  has  proved  that  under  appropriate 
conditions  a  consistent  scheme  is  convergent  if  and  only  if  it  is  stable. 

The  Lax  Equivalence  Theorem  states  that  given  a  properly  posed  linear 
initial  value  problem  and  a  finite  difference  approximation  to  it  that  satisfies  the  con¬ 
sistency  condition,  stability  is  the  necessary  and  sufficient  condition  for  convergence 
(see  [Ref.  15]). 


4.  Further  Examples 

Given  a  PDE  and  a  finite  difference  mesh  one  can  use  any  of  the  following 
procedures  to  develop  a  finite  difference  scheme. 

a.  Tables 

b.  Taylor  series  expansions 

c.  polynomial  fitting 

d.  integral  methods 

e.  control  volume  techniques. 

One  may  get  the  same  scheme  by  using  different  approaches.  As  an  example 

^  r 

for  procedure  h  we  develop  a  three  point  second  order  approximation  for  —  on  a 
df 

nonuniform  mesh.  —  at  point  0  can  be  written  as  a  linear  combination  of  values  of 
/  at  A,  O,  and  B, 


dx 


=  CJ{A)  +  C2f{0)  +  C3f{B) 


Figure  1.  Nonuniform  mesh 

We  use  Taylor  series  to  expand  f{A)  and  f{B)  about  the  point  O, 
m  =  no  -  Ax)  =  /(O)  -  Axf'iO)  +  ^/"(O)  -  ^/'"(O)  ± 

2  D 
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/(B)  =  f(0  +  aAx)  =  /(O)  +  aAx/'(0)  +  ^|-/''(0)  +  —^f"'(0)  +  ■  •  • 
Thus 


dx  o 


{C\ C2  Cz)f{0)  +  (aCs  —  Ci)Aa;— I  +  {Ci  +  cx^Cs)  ^ 

+  (a"C3-C.)^^|  +■■■ 


This  yields  the  following  system  of  equations 


Cl  +  C2  +  C3  —  0 

-C.  +«C3  =  ^ 

Cl  +^2(73  =  0 

The  solution  is 

^  _ _ ^ ^  ^  ^  ^ _ 

^  (a  +  l)Aa:’  ^  ccAoj  ’  ^  a(Q:  +  l)Aa: 

and  thus 

df  -off (A)  +  (o^  - 1)/(0)  +  /(B)  a 

5a:  a(Q;  +  l)Ax  6  dx^  o 

Note  that  if  the  grid  is  uniform  then  a  =  1  and  this  becomes  the  familiar  centered 

difference. 


5.  Irregular  Mesh 

It  is  more  convenient  to  use  a  uniform  mesh  and  it  is  more  accurate  in  some 
cases.  However,  in  many  cases  this  is  not  possible  due  to  boundaries  which  do  not 
coincide  with  the  mesh.  In  this  case  several  possible  cures  are  given  in  [Ref.  14].  The 
most  accurate  of  these  is  a  development  of  a  finite  difference  approximation  which  is 
valid  even  when  the  mesh  is  nonuniform.  It  can  be  shown  that 

^  2  /Ue  -  up  _  Up  -  ua 

p  [1  4”  ^  Gch^  /lx 

Similar  formula  for  Uyy.  Note  that  for  0=1  one  obtains  the  centered  difference 
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Figure  2.  Irregular  mesh  near  curved  boundary 


approximation. 

Due  to  the  need  to  refine  the  mesh  in  some  of  the  domain  to  maintain  the 
accuracy  one  is  advised  to  use  a  coordinate  transformation  which  is  covered  in  Neta 
(see  [Ref.  13]). 


6.  Stability 

The  problem  of  stability  is  very  important  in  numerical  analysis.  There  are 
two  methods  for  checking  the  stability  of  linear  difference  equations.  The  first  one 
is  referred  to  as  Fourier  or  von  Neumann  and  assumes  the  boundary  conditions  are 
periodic.  The  second  one  is  called  the  matrix  method  and  takes  care  of  contributions 
to  the  error  from  the  boundary. 

a.  von  Neumann  analysis 
Suppose  we  solve  the  advection  equation 


dt  dx 


by  Lax  method  (see  [Ref.  13]) 


1 


At 
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If  a  term  (a  single  term  of  Fourier  and  thus  the  linearity  assumption) 


F 

1  tj 


_  gCit gZ/u77iir 


is  substituted  into  the  difference  equation,  one  obtains  the  amplification  factor 

=  cos  13  —  iu  sin  (3 


where 


Courant  number 


(3  =  km^X  . 


The  stability  requirement  is 


e“‘|  <  1 


and  implies 


M<i. 


This  is  called  the  Courant-Friedrichs-Lewy  (CFL)  condition. 

b.  Matrix  method 

Suppose  again  we  solve  the  advection  equation  using  Lax  method  but  now  we  assume 
periodic  boundary  conditions,  i.  e. 


m+ln 


=  F, 


In 


The  system  of  equations  obtained  is 


£n+l  = 


where 
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A  = 


0 

^0 

0 

0  0 

1-1/ 


1  -i/ 


0 

1-1/ 


0 

1+iy 


l  +  iy  1 


0 

1  -i/ 


It  is  clear  that  the  eigenvalues  of  A  are 


\  Stt  ^ ,  Stt  . ,  ,  , 

Aj  =  cos — (J-l)^-^l/sln — U-1),  ^  =  1,' 

m  m 

Clearly  the  stability  of  the  method  depends  on 


,  m . 


\P{A)\<1. 


Note  that  one  obtains  the  same  condition  in  this  case.  The  two  methods  yield  identical 
results  for  periodic  boundary  condition.  It  can  be  shown  that  this  is  not  the  case  in 
general.  See  work  by  Hirt  (see  [Ref.  16]),  Warming  and  Hyett  (see  [Ref.  17])  and 
Richtmeyer  and  Morton  (see  [Ref.  15]).  We  will  explore  the  stability  of  the  shallow 
water  equations  in  more  detail  using  Fourier  techniques  in  the  next  chapter. 

7.  Example:  Shallow  Water  Equations  in  2D 

Arakawa  and  Lamb  ([Ref.  18])  have  investigated  a  finite  difference  scheme  for 
the  nonlinear  shallow  water  equations  using  square  and  staggered  square  grids.  We 
will  show  two  of  their  examples,  one  of  a  square  unstaggered  grid,  and  one  example 
of  a  staggered  square  grid.  For  more  information  on  finite  difference  schemes,  refer 
to  their  work.  [Ref.  18] 
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a.  Square  Grid 

Recalling  the  linearized  version  of  the  shallow  water  equations 
du 


dv  dh  ^ 


dh  tj  ( 

dt  ^  \dx^  dyj 


=  0 


? 


Arakawa  and  Lamb  (see  [Ref.  18])  use  the  following  finite  difference  approximation 
for  the  unstaggered  square  grid. 


^  +  fn  +  ^(S,hy  =  0 

^ + 7  {(u^r  +  (^r)  =  0 


where  {SxliY  =  in  our  notation  this  is  fixSxh. 

Ci 

b.  Staggered  Square  Grid 

In  siniulation,  according  to  Arakawa  and  Lamb  (see  [Ref.  18]),  a  stag¬ 
gered  grid  approach  is  best.  The  reason  that  a  staggered  square  grid  is  better  for 
modelling  shallow  water  flow  than  an  unstaggered  square  grid  is  that  for  the  same 
time  step,  tlie  staggered  grid  will  give  higher  levels  of  accuracy.  Now  let  us  look  at 
Arakawa  and  Lamb's  example  of  a  staggered  square  grid.  Recalling  the  shallow  water 
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equations  (2.10). 


du  du  du 

dh 

= 

dv  dv  dv 

dh 

dh  d  ,  ,,  d  , 

■57  +  5; 

=  0. 

(4.1) 


Notice  that  we  substituted  h  for  {h  —  hs)  therefore  ignoring  the  bottom  topography. 
For  the  staggered  grid,  following  Arakawa  and  Lamb’s  example  (see  [Ref.  18]),  multi¬ 
ply  the  first  of  (4.1)  by  h  and  add  it  to  the  last  of  (4.1)  multiplied  by  u,  also  multiply 
the  second  of  (4.1)  by  h  and  add  it  to  the  last  of  (4.1)  multiplied  by  v.  This  gives  us 

hut  +  huux  +  hvuy  —  //iu  -f  ghhx  4-  uht  4-  u{hu)x  +  u{hv)y  =  0 
hvt  4-  huvx  4-  hvvy  +  fhu  4-  ghhy  4-  vht  4-  v{hu)x  4-  'v{hv)y  —  0. 

This  reduces  to 

{uh)t  4-  {huu)x  4-  {hvu)y  —  fhv  +  ghh^  =  0 
{vh)t  4-  {huv)x  4-  {hvv)y  4-  fhu  4-  ghhy  =  0. 

This  is  another  useful  form  of  the  momentum  equations.  Now  multiply  the  first  of 
(4.1)  by  uh  and  add  it  to  the  last  of  (4.1)  multiplied  by  ,  also  multiply  the  second 

Li 

of  (4.1)  by  vh  and  add  it  to  the  last  of  (4.1)  multiplied  by  — .  This  gives  us 

IZ^ 

uhut  4-  hu^Uj  4-  huvuy  —  fhuv  +  ghuh^  4-  —ht  4-  -r{hu)x  4-  -7r{hv)y  =  0 


vhv,  4-  hvuvx  4-  hv^Vy  4-  fhvu  4-  ghvhy  4-  y/^t  4-  y (^«)x  4-  y (/ii’)y  =  0. 
This  now  reduces  to 

(^^y)  ~  +  dMx  =  0 


fiy^  4-  f/iuy^  =  0- 


(4.2) 
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This  gives  the  equation  for  the  time  change  of  kinetic  energy  according  to  Arakawa 
and  Lamb  (see  [Ref.  18]).  If  we  multiply  the  last  of  (4.1)  by  gh,  we  get 

gh  [ht  +  {hu)x  +  {hv)y]  =  0 


or 

fl'y)  +  {gh'^u)x  +  {gh^v)y  -  gh{uhx  +  vhy)  =  0.  (4.3) 

Note  that  the  Coriolis  force  does  not  contribute  to  the  change  in  total  kinetic  energy. 
.Also  note  that  the  sum  of  terms  in  (4.2)  and  (4.3)  is  zero.  It  then  follows  that  we 
have  conservation  of  energy  which  can  then  be  used  in  the  construction  of  the  finite 
difference  scheme. 

Arakawa  and  Lamb’s  choice  for  the  differencing  of  the  continuity  equa¬ 
tion  (see  [Ref.  18])  was  decided  in  an  effort  to  keep  it  as  simple  as  possible.  At  h 
points,  the  last  equation  of  (4.1)  can  be  represented  as 


wliere 

=  d[hyu],^ij 

are  the  mass  fluxes  and  are  defined  at  u  and  v  points,  respectively.  This  is  semidiscrete 
therefore  the  time  change  terms  will  be  left  in  differential  form. 

For  the  finite  difference  scheme,  the  total  kinetic  energy  is  conserved 
during  the  inertial  process.  Therefore  the  terms 


{uh)t  -f  {huu)x  -f  {him)y 


ran  In*  written  as  follows 

a  .....  .  1 


We  will  use  Arakawa  and  Lamb’s  reasoning  to  define  and 

in  the  next  few  paragraphs  (see  [Ref.  18]).  (i,j)  are  used  as  the  indices  and  are 
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chosen  in  such  a  way  that  they  satisfy 


+  4'^^“^] . .  =  0 

Now  multiply  (4.5)  by  Uij  and  subtract  it  from  (4.4),  we  will  get 


(4.5) 


r(M) 


1  h 


H\y^  +  jp  +  ^1(“)4«2'  +  .7"W4u^'  +  G(^)5'yuy' 


(4.6) 


1 


Take  (4.6)  and  multiply  it  by  Uij  and  add  it  to  (4.5)  which  is  multipied  by  ^ufj  and  a 
finite  difference  analog  of  the  first  three  terms  of  the  first  equation  of  (4.2)  is  obtained. 
The  finite  difference  analog  is  written  the  same  as  Arakawa  and  Lamb  wrote  it  and 
is  as  follows  (see  [Ref.  18]). 


1 


{u) 


2(PV 


I  nM 

J+1 

2  2  2  2 

Notice  that  in  this  equation,  the  kinetic  energy  flux  term  appears  twice,  but  with  the 
opposite  sign  the  second  time.  The  definition  of  these  terms  is  not  dependent  on  the 
form  of  the  equation  as  long  as  the  total  kinetic  energy  does  not  increase  or  decrease 
over  the  domain. 

The  Coriolis  term  —fliv  can  be  represented  at  u,j  by 


and  fhu  at  by 
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The  pressure  gradient  terras  ghh^  and  ghhy  can  be  represented  as  g[h^5xh]ij  and 
g[h'^Syh]-_^i.j^i  respectively.  This  is  a  general  momentum  advection  scheme  for  non- 
divergent  flow  to  a  scheme  that  will  maintain  conservation  of  total  energy  and  also 
maintain  the  divergent  flow.  Horizontal  errors  due  to  discretization  should  be  small 
on  the  planetary  scale  due  to  the  large  size  of  the  system  relative  to  the  grid  size. 
For  a  more  detailed  rendition  of  this  example,  refer  to  Arakawa  and  Lamb  (see  [Ref. 
18]).  Neta  and  Lustman  generalized  this  to  nonuniform  grid  (see  [Ref.  19])  . 

B.  FINITE  ELEMENT 

In  the  finite  element  methods  (FEM),  the  domain  is  divided  into  subregions 
railed  finite  elements,  hence  the  name.  The  unknown  function  u  is  represented  as  the 
interpolating  polynomial  in  each  element.  This  representation  is  continuous  along 
with  its  derivatives  (to  a  certain  order)  in  each  elment. 

1.  Basic  Concepts 

One  of  the  ways  to  formulate  the  finite  element  is  via  the  so  called  weighted 
residuals  method.  In  this  method,  the  desired  function  u  is  replaced  by  a  finite  series 

N 

u’^  = 

i=i 

The  set  of  functions  are  called  basis  functions.  Clearly  one  can.  not  expect  to 
satisfy  the  partial  differential  equation, 

Lu  =  f. 

The  re.sidual  R  is  defined  as 

R  =  Ln'*  -  /. 

In  order  to  obtain  the  undetermined  coefficients  Uj,  one  sets  the  weighted  residuals 
to  zero,  i.e. 

J  Rwi  =  0 
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where  the  weights  to,-  may  be  chosen  in  various  ways.  If  the  weights  are  chosen  to  be 
the  same  as  the  basis  functions,  one  obtains  Galerkin’s  method.  For  collocation,  the 
weights  are  Dirac  delta  functions.  In  the  subdomain  method  one  uses  the  character¬ 
istic  function  of  each  subdomain  as  the  weight,  i.e.  tOj  =  1  in  the  subdomain  D,-  and 
zero  elsewhere. 


Example 

Given  the  equation 

(Pu 

-r—  +  u-\-x  =  0  0<x<l 

dx^ 

with  homogeneous  boundary  conditions 

u(0)  =  'u(l)  =  0, 

we  suppose  that  we  can  approximate  u{x)  by 

=  ai^i{x)  -I-  a2<l>2{x) 

and  let  the  basis  functions  be 

(f>i  =  x[\  —  x) 

(4.7) 

^2  =  —  x). 

Notice  that  each  basis  function  satisfies  the  boundary  conditions  and  therefore 
satisfies  those  also.  The  residual  is 

R  =  ~YT  +  a:  =  (—2  -f-  x  —  x^)ai  +  {2  —  6x  +  x"^  —  x^)a2  +  x.  (4.8) 

Tlie  collocation  method  yields  (using  x^  =  ^,X2  =  \  as  collocation  points)  the 
following  two  ecpiations  for  the  unknowns  aj  and  02. 

_  35  _  1 

ifiOi  g^a2  —  ^ 

4®!  +  =1 
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The  approximate  solution  is  then 


u 


x{\  —  x) 


42  +  40a: 
217  ■ 


The  subdomain  method  yields  the  following  system  when  the  interval  is  sub¬ 
divided  into  two.  ^ 

Rdx  =  0 

Jo 


/; 


Rdx 


=  0 


11 


53 


'12“^  192 


02 


1 

'8 


11 


229 


- U] - 

12  192 


02  =  -- 


(4.9) 


Ol 


97 

517 


02 


88 

517 


u 


=  a:(l  —  a:) 


97  +  88x 
517 


For  Galerkin’s  method  we  have  the  following  two  integrals  to  evaluate. 


/  R(f)i  dx  = 
Jo 

R(f>2dx  = 
0 


(4.10) 


with 


3  3 

To“'  +  10“" 


1 

12 


3  13  1 

70°' Tos"’  “  20 
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where 


71 

369 


02 


7 

41 


b,  /I  \/  \ 

“  =  +  41^) 

The  accuracy  of  the  approximate  solution  depends  on  the  choice  of  the  basis 
functions  and  the  method  used. 

2.  Weak  Formulation 

The  fundamental  integral  statement  of  the  finite  element  methods  can  be  in¬ 
terpreted  as  a  combination  of  a  weighted  residual  and  a  process  of  integration  that 
reduces  the  order  of  continuity  required.  Going  back  to  the  previous  example 
<Pu 


u(0)  =  u(l)  =  0 


+  a2<i>2{x) 


(4.11) 


/■ 


<Pu 


[dx'^ 

Integration  by  parts  yields 
du  d4>i 


X 


fy 


(j)idx  =  0 


du 


*  =  1,2 


+  u(j),  +  xcf>i)dx  +  —(j)i 
ax  ax  ax 


=  0 


1,2 


It  is  in  this  weak  form  that  one  substitutes  the  approximate  solution.  Notice  that 
<p,(0)  =  0,(1 )  =  0  and  thus  the  boundary  term  vanishes. 


3.  Choice  of  Basis  Functions 

The  accuracy  of  the  methods  of  weighted  residuals  depends  mainly  on  the 
choice  of  basis  functions.  In  the  previous  example,  we  have  chosen  polynomials.  We 
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now  introduce  several  possible  piecewise  polynomials.  The  simplest  polynomial  is 
piecewise  constants. 

a.  Piecewise  Constant 


b. 


4)i{x) 


I  a  xe(xi-i,xi+i) 

0  elsewhere 


Piecewise  Linear 


f  X-Xi-i 


Xi 


Xi—\ 


Xi-1  <  X  <  Xi 


(j)i{x)  =  < 


Xj+i  -  X 
Xi 


Xi  <  X  <  ajj+i 


[  0  elsewhere 

c.  Piecewise  Quadratic 

For  higher  order,  it  is  easier  to  introduce  a  nondimensional  coordinate 
^  where  (  —  1  <  ^  <  1).  This  choice  facilitates  numerical  integration  by  Gaussian 
<|uadratures. 


■fS-.K)  =  -ie(i-f) 


MO  =  i-e 
MO  = 

The  subscript  indicates  the  point  (  at  which  (/>  is  one.  o  is  zero  at  the 
other  two  points  as  shown  in  the  following  figure. 


38 


Figure  3.  Piecewise  quadratic  basis  functions 
d.  Piecewise  Cubic 

It  is  easy  to  see  that  piecewise  cubic  basis  functions  are  given  in  the 
nondiinensional  coordinate  ^  by 

4>-,(0  =  i(i-f)(9f^-i) 

- 1)  ■ 

MO  =  -^(1  +  e)(9f  - 1). 

The  above  basis  functions  interpolates  along  the  element  using  a  Lagrange  third  order 
polynomial.  Another  choice  is  the  Hermite  polynomial.  In  addition  to  continuity 
of  second  derivatives  over  the  element,  one  has  first  derivative  continuity  between 
elements.  Therefore  Hermite  polynomials  interpolate  the  derivatives  also.  This  is 
useful,  for  example,  in  flow  problems  where  one  has  to  differentiate  the  potential  to 
oiitain  the  velocity  field. 
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Each  node  is  identified  with  two  basis  functions. 

■A>.  =  |(i-0"W  +  2)  <A„(-i)  =  i,  *,,(!)  =  .^;,(±i)  =  o 

4*02  =  ^2(1)  =  1,^2(  —  1)  =  *^02(^0  “  ^ 

<^'.1  = 

(*'12  =  + 


4.  2D  Basis  Functions 

It  is  very  easy  to  extend  the  method  of  weighted  residuals  to  higher  dimensions. 
In  this  section  we  discuss  several  possibilities  of  bcisis  functions  for  rectangular  and 
triangular  elements.  We  close  by  discussing  isoparametric  finite  elements  (irregular 
quadrilateral  elements). 

a.  Rectangles 

As  in  one  dimension,  one  can  present  Lagrangian  and  Hermitian  basis 
functions.  The  Lagrangian  basis  functions  are  obtained  when  a  product  of  two  one 
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dimensional  basis  functions  are  used.  We  will  use  a  quadratic  for  this  example. 


<f>l  = 

-0(1  -^) 

(-1.-1) 

<f>2  — 

(0,-1) 

<t>3  = 

-^^(1  +  0(1  -0 

(l.-l) 

<l>4  — 

-101-0(1-^0 

(-1.0) 

(t>h  = 

(1  -  0)(1  -  ^0 

(0,0) 

<f>e  = 

101  +  0(1-^0 

(1,0) 

<f>7  = 

-KOI  -0(1  +  0 

(-1,1) 

<f>s  = 

I(l  -0)»7(l  +  0 

(0,1) 

<f>9  = 

K'i(i  +0(1+0 

(1,1) 

Figure  5.  Nodes  associated  with  a  rectangular  element 
The  9  nodes  associated  with  the  rectangular  element  are  shown  in  the 
figure  above.  0,  is  the  basis  function  having  a  value  of  one  at  the  point  listed  next  to 
its  definition  and  zero  at  the  other  8  nodes. 
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It  is  easier  to  write  the  basis  functions  by  distinguishing  between  corner 
nodes,  side  nodes,  and  interior  nodes. 

Quadratic: 


Corner  node 

+  ^6)Wi(l  +VVi) 

Side  node. 

6  =  0  lTiT]i{l  +  r]r]i){l  - 

Side  node. 

^t  =  0  1^6(1 +^6)(i 

Interior  node  (1—  —  77^) 

Linear: 

^(1  H-  +  CCi) 

Another  possibility,  which  preserves  the  properties  of  the  Lagrangian 
basis  functions,  but  eliminates  the  interior  nodes  are  called  serendipity  basis  functions. 
For  example,  in  the  case  of  the  quadratic: 

Corner  node  i(l  +  +  ■nr]i)i^(i  +  VVi  “  1) 

Side  node,  =  0  |(1  -  ^2)(1  +  ■qTji) 

Side  node,  77i  =  0  |(1  -  ?72)(1  + 

The  cubic  and  Hermite  cubic  can  be  found  for  example  in  Lapidus  and 
Finder  [Ref.  20]. 

6.  Triangles 

The  triangular  element  is  the  most  well  known.  It  allows  more  accurate 
presentation  of  an  irregular  domain  than  the  rectangular  element.  A  natural  repre¬ 
sentation.  calh'd  area  coordinates  system  was  introduced  to  simplify  the  integration. 


Figure  6.  Area  coordinates 


Each  coordinate  Li  is  defined  in  terms  of  the  triangle  sub  area  whose 
base  is  the  line  L,-  =  0  and  the  vertex  is  the  point  /,•  =  1,  thus 


where  A,  are  the  areas  of  the  subtriangles  and  A  is  the  area  of  the  whole  element. 
C'learly  A  =  Ai  +  A2  +  A3.  Also  clear  is  the  Li  is  unity  at  node  i  and  zero  at  the 
other  two  nodes.  For  example  =  1  if  P(x,y)  coincides  with  the  vertex  (xi,j/i)  and 
it  is  zero  if  P  coincides  with  either  of  the  other  two  vertices. 

Note: 


Li  L2  L3  —  1 


milmaima! 


/  L'T'L'^^L^^dxdy  =  2A- 
Jt  ^  ^  ^  ^  {Tm+m2  +  m3  +  2)\ 

where  77? m2,  m3  are  nonnegative  integers. 

It  dy  dy^^^^  ~  ~  ^t+i)(3Jj+2  —  ^j+i) 


I 


T  Ox  dx  ~  ~  J/*+2)(yj+i  -  yj+2) 


where  ;  and  j  are  ryclic  for  1  <  i.j  <  3. 

In  many  cases  one  requires  or  prefers  numerical  integration.  We  will 
not  disruss  this  matter,  but  the  reader  is  refered  to  Connor  and  Will  (see  [Ref.  21]), 
Br<‘bbia  and  Connor  (see  [Ref.  22])  et  al. 


For  quadratic  basis  functions,  one  requires  three  more  nodes  per  ele¬ 
ment  and  these  are  taken  cis  midpoints  of  each  side.  The  six  basis  functions  in  area 
coordinates  are 

4Ll  - 
ALi  L2 
2Ll  -  L2 
4L2L3 
2Ll  -  Ls 
4L3T1. 

Several  possible  higher  order  triangular  elements  can  be  found  in  Fe- 
lippa  [Ref.  23],  Lapidus  and  Finder  [Ref.  20],  and  Zienkiewicz  [Ref.  24]  and  others. 

c.  Isoparametric 

In  the  Isoparametric  finite  elements,  all  irregularly  shaped  elements  are 
mapped  onto  regular  elements  in  order  to  help  with  the  integration. 

Examples 


1.  Any  quadrilateral  with  straight  sides  can  be  mapped  onto  a  square 
whose  vertices  are  at  (±1,±1). 

2.  Any  quadrilateral  with  curved  sides  can  be  mapped  onto  the  same 

square. 

3.  Any  quadratic  triangle  with  curved  sides  can  be  mapped  onto  an 
equilateral  triangle. 


It  can  be  shown  the  the  mapping  from  the  irregular  shape  in  x^y  coor¬ 
dinate  to  (^,7/  is  given  by 


=  E.=,:r,-0.(e,7/) 
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where  {xi^yi)  are  the  vertices  of  the  quadrilateral  and  <l>i{^,r])  are  the  basis  functions 
of  the  appropriate  order,  (in  the  first  case  the  <f>i  are  linear.) 

The  evaluation  of  integrals  require  the  Jacobian  of  transformation  J. 

'  dx  dy  ' 

J  = 

dx  dy 
.  dy  dy  . 

Isoparametric  elements  of  higher  order  can  be  found  in  Lapidus  and 
Finder  (Ref.  20],  Zienkiewicz  [Ref.  24]  and  others. 

5.  3D  Basis  Functions 

The  extension  to  three  dimensions  is  straight  forward  if  one  has  hexahedral 
elements  either  Lagrangian  or  serendidpity  (no  nodes  in  the  interior).  One  can  also 
use  tetrahedral  or  pentahedral  elements. 


Figure  7.  Three  Dimensional  Elements 

Tables  for  serendipity  type  can  be  found  in  Zienkiewicz  [Ref.  21].  \’erge  [Ref. 
25]  et  al.  Area  coordinates  are  replaced  by  volume  coordinates. 
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6.  Example:  Shallow  Water  Equations  in  2D 

In  this  section,  we  give  an  example  of  solving  the  shallow  water  equation  using 
finite  elements.  We  will  ignore  the  bottom  topography.  The  shallow  water  equations 
in  vorticity  divergence  formulation  are  (as  one  may  recall  from  Chapter  2) 

C<  +  ('^Q)x  +  ('^Q)y  =  0 

Dt  +  gV%  +  V^K-(vQ)x  +  (uQ)y  =  0 

fit  +  [uh]^  +  [vh]y  =  0, 

where  h  is  the  geopotential  height,  u  the  east/west  component  of  the  wind,  v  the 
north/south  component  of  the  wind  and  /  is  still  the  Coriolis  parameter.  Kinsman 
(see  [Ref.  26])  shows  that  the  equations  can  be  discretized  using  finite  elements  as 
follows: 


.  dVj  dVi  ,  .  dVj  dVi  ,  \  f\d  , 


dx  dx 


dy  dy 


dx 


Vi 


/ 


gy  k  gy 


Vi  + 


,  dVi  ,  , ,  ,  avA ,, 


/  + /  it  h  - 


N 

s 


dVj  dVj 

dt  dx  dx  dt  dy  dy 


(4.12) 


(4.13) 

_  [  ^ 

J  \  dt  dx  dx  dt  dy  dy  dt  dx  dx  dt  dy  dy  j 

Jtx  -  /I;  ((“0)'^  - 

(4.14) 


46 


du  dv 


dcf) 


Divergence  (V^x)  =  ^  ^  'ir  “  “'“/o  along  the  boundary  walls.  For  further 

ox  oy  oy 

look  at  this  example  see  [Ref.  26].  We  have  used  Vi  for  the  basis  functions  in  order 


not  to  confuse  it  with  the  geopotential  height.  Note  that  (4.13)  and  (4.14)  are  time 
dependent.  Also  note  that  Kinsman  is  using  a  leap-frog  time  discretization  to  obtain 


each  next  time  step.  Leap-frog  requires  an  additional  starting  value  obtained  by  the 
Matsuno  scheme.  (See  [Ref.  5].) 
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V. 


STABILITY  ANALYSIS 


Before  getting  into  the  stability  analysis  of  the  shallow  water  equations,  this 
paper  shall  review  some  basics  about  integral  transforms  by  focusing  on  the  Fourier 
transforms.  Integral  transforms  are  applied  in  mathematical  type  operations  by  tak¬ 
ing  an  equation  which  is  unsolvable  (or  really  hard  to  solve)  in  its  original  form  and 
transforming  it  to  an  equation  that  is  algebraically  solvable.  Generally,  the  harder 
the  original  problem  is  to  solve,  the  harder  the  inverse  transform  is  to  find.  This 
chapter  will  not  go  into  the  nitty-gritty  of  integral  transforms,  but  will  offer  some 
insight  into  why  and  how  they  work  before  applying  transforms  to  the  problem  at 
hand.  For  more  information  on  integral  transforms  see  Miles  [Ref.  27].  Simply  put, 
an  integral  transform  is  of  the  form  J^{p)  =  f^  K{p,x)f{x)dx  where  the  function  to 
be  transformed  is  f{x),  and  K{p,x)  is  the  transforming  function  which  is  known  ais 
the  kernel  of  the  transform. 

A.  ONE  DIMENSION 

1.  Fourier  Transform 

The  Fourier  Transform  in  one  dimension  is 

f{k)  =  T{f{x))  =  j  e-^’^^f{x)dx. 

J  — CX5 

Also,  recall  that  the  inverse  Fourier  Transform  is 

/(x)  =  r-'(!{k))  = 

2.  Shallow  Water  Equations 

Srhoenstadt  (s<*c  [Ref.  2S])  studied  the  effect  of  replacing  the  spatial  deriva¬ 
tives  in  a  dispersive  wave  equation  by  using  Fourier  transform  techniques.  We  will 
repeat  some  of  his  work  here,  but  for  a  more  in  depth  understanding,  the  reader 
should  consult  the  original  text.  Let  us  consider  the  effect  of  semi-discretization  on 
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the  solution  of  a  shallow  water  type  model.  The  shallow  water  equations  in  a  one¬ 
dimensional  linearized  form  with  no  mean  flow,  as  one  may  recall  from  Chapter  3, 
are 

ut-  fv  +  ghx  =  0 

vt  +  fu  =0  (5-1) 

hfi  -j-  fIxLx  —  0. 

Now  if  we  use  the  Fourier  transform  on  these  equations,  (hats  ( " )  symbolize  Fourier 
transforms),  we  arrive  at 

Ut  =  fv  —  ikgh 
vt  =  -fu 
hi  =  —ikHu 

The  initial  conditions  are 

uo  =  u(k,0)  =  «(a:,0)e”^*®dx 

Vo  =  ^(A;,  0)  =  f^v(x,0)e~‘^'^dx 

ho  =  h(k,0)  =  fZo  /i(a:,0)e-‘*^dx 

The  transformed  shallow  water  equations  are  a  coupled  set  of  constant  coefficient 
ordinary  differential  equations  which  can  be  solved  fairly  easily  (see  [Ref.  28])  and 
are 

-  iut  ^3  -iut 

u  =  — e - e 

u  V 


V  =  ikgai  -t- 


*7"2  iut  ,  *7o;3  -l 


€'''■  + 


tvi 


(5.2) 


h  =  /oi  - 


kHa2  .  kHas  ■ 


d- 


xi/i 


whrrr 


-2  _ 

f  ■ 


A  = 


1  he  o,’s  are  picked  to  satisfy  the  initial  conditions.  Recalling  that  sinx  = 


2i 


^IT  ^  g-ll 

aiid  COST  = - we  can  rewrite  (5.2)  by  using  the  intial  conditions  to  solve  for 
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O',  and  then  simplifying  to  get 


w,  ,x  .  .  ,  f^o  .  .  ikgho  . 

u[K,t)  =  uq  COS  ut-\ - smut - smut 

u  u 


v[k,t)  = 


fuo  .  ^  ,  f  k^Hf  p  A  n  i 

sin  ut  +  <  — - 1 — -  cos  ut>  Vo  -I - —  {1  —  cos  ut}  ho 


u 


ikHuo  .  ikHf 


P  .  UHf 


h{k^  t)  = - sin  i/t - —  {1  —  cos  ut}  vo  +  <  H - cos  ut  >  ^0 


u 


u^ 


u^  u^ 


The  steady  state  solutions  are 


or  to  rewrite 


Us{k)  =  0 


PgH ,  ikgf  f 
=  -^^o  +  -^ho 

t  (U\  -  f't, 

hs{k)  —  ^2  ^0  ^2^° 


Us{k)  =  0 


-  (h\  -  A.  ^ 

Vs{k)  —  Uo  H - r  I  —^ho  —  Vq 

\  f  } 


By  noting  that 


fts(fc)  =  Ao  +  (^^^0  -  yoi 


/. 


00  _ 


-^  ikxd  _  - ^  ‘^AIl 

00  1  +  A:2A2  1.2 
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and  using  the  convolution  theorem,  Schoenstadt  (see  [Ref.  28])  inverted  the  steady 
state  transformed  equations  to  yield 

Us{x)  =  0 

v,(a;)  =  v{x,  0)  +  ^  (fix 

hs{x)  =  h(x,  0)  -  2^  sgn(x  -  s)e-^  (f  H 

Note  that  in  the  one  dimensional  case,  u{k,  0)  does  not  contribute  to  the  steady  state 
solutions.  Schoenstadt  (see  [Ref.  28])  goes  on  to  explore  the  transient  parts  of  the 
model,  and  these  can  be  studied  in  more  detail  in  his  1977  paper. 

Schoenstadt  (see  [Ref.  29] )  also  analyzed  a  difference  scheme  for  one  dimension 
shallow  water  equations  that  is  similar  to  the  one  we  will  present  later  in  this  chapter 
for  two  dimensions. 

B.  TWO  DIMENSIONS 
1.  Fourier  Transform 

Recall  that  the  Fourier  Transform  in  2  dimensions  is 

hv,„w,)  =  :F(f{x,y))  =  r  r 

J—oo  J—oo 

and  the  inverse  Fourier  Transform  is 
/(X,y)  = 

Also  recall  several  elementary  properties  such  as: 

^(/.)  = 

^ifx)  = 

J^ify)  =  iw2J^{f) 

=  iwTif) 

^(vV)  = 
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We  will  use  several  of  these  properties  later  in  this  chapter  when  dealing  with  the 
continuous  case  and  also  the  semi-discrete  case.  The  convolution  theorem  will  also 
be  used  in  this  section  and,  in  two  dimensions,  it  is 

^  poo  roo 

f(x,y)  =  J^~'^{g{k,l)-h{k.l))=  /  g{xo,yo)h(x  -  xo,y  -  yo)dxodyo. 

j-ooj’-oo 

2.  Continuous  Case 

Recall  the  linearized  shallow  water  equations  from  Chapter  3. 

du 


dv  dh 


(5.3) 


ah  „  fBu  dv\  „ 

Following  the  work  of  Neta  (see  [Ref.  2]),  we  will  go  through  the  reasoning 
done  by  Neta  and  by  Neta  and  DeVito  (see  [Ref.  7])  to  arrive  at  their  solution.  Let 
us  take  the  first  equation  of  (5.1)  and  transform  it  using  Fourier, 


( du  dh\  [°°  ( du  .  dh\  ,  , 

(a  - j  =  LL  ( aF  ■  j " 


=  ^,r  r  -fTr  +  gr  r 

Ut  J—ooJ  ^oo  J—ooJ—oo  J—ooJ—ooUX 


dh 


du 

=  —  -  fv  -f  igkh  =  0. 


Similarlv 


[dt 


+  9^  ^9lh  =  0. 


\dt  '  "  \dx  '  dy))  dt 


=  TT  +  iHku  4-  iHlv  =  0. 
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variables  used  there.  They  are 


D 

c 

dx 

d<,i 


=  ikuo  +  iivo 

=  Uuq  —  ikvo 

=  ijkho  -  Vo 
=  iylho  +  Uo 


where  the  steady  state  solution  is 


Us 

Vs 


—  Uo  -{■  j^{X^ikD  —  dy) 

=  H — ~{X^ilD  +  dx) 


H  p  .  ^ 

Jxg  —  ho  “h  ~~z  ~(tkdx  "h  zldy^ 

f 

We  can  now  invert  these  equations  using  the  convolution  theorem  and  the  following 
integrals. 

1  /‘27r 

—  /  =  Jo(rp) 

ZTT  JO 


1  /•«>  Jo(r/9) 


27r 


too 

Jo  T 


27rA2 


(5.6) 


-^(t\  = 

wJ 


1  f 


27rA2 


J()  and  A*o  are  Bessel  functions  of  order  zero.  Using  these  integrals,  the  inverse 
obtained  from  them,  and  the  convolution  theorem,  we  have  the  steady  state  solutions. 
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u 


\  too  too  Q 

..{x,y)  =  ‘^o{x,y)  - —J^^J  —Ko(q{(T,T))D{(7,T,0)dadT 

too  too 

■/  /  KQ{q{cr,T))dy(a,T,0)dadTj 

J—ooJ—oo 


27rA2 


2  ACO  fOO  Q 

Vs{x,y)  =  V(i{x,y)  -  —  J^^J^^~Koiq{(T,T))D{ar,T,0)dadT 

too  too 

-/  /  Ko{q{(T,T))d::,{(T,T,0)d(TdT, 

'  t/— ooJ  — oo 


and 

hs{x,y)  =  ho{x,y)- 


2ttX^  j  — oo»y  — oo 

H  1 


/27rA2 


noo  /  /) 


d 


+—Ko{q{(7,r))da:{a,T,0)j  dadr 


where 


q{a,T)  = 

^Jix-af  +  {y-Tf 

A 

D{a,  T,  0)  = 

duo  dvo 
da  ^  dr. 

4(a,T,0)  = 

gdho 
f  dr 

a 

3 

II 

gdho  , 

3.  Semi-discrete  Equations 

Here  we  will  present  the  semi-discrete  shallow  water  system  and  then  show  the 
inverse  Fourier  transform  for  certain  selections  of  q,  /?,  and  7’s.  Semi-discrete  means 
that  time  is  not  discretized.  In  the  semi-discrete  case,  the  shallow  water  equations 
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are 


-  /3fv  +  ig'fih  =  0 


+  f^fu  +  ig-f2h 


0 


(5.7) 


a—  +  iH  (7iu  +  72v)  =  0 

where  a,  /i,  71,  and  72  depend  on  the  scheme  used  (see  tables  II  and  III)  from  Neta 
[Ref.  30]. 


Table  II.  Filter  weights  for  second  and  fourth  order  finite  differences 


Scheme 

a 

7i 

72 

A2 

1 

1 

sin  kd 

sin  Id 

d 

d 

B2 

1 

1 

sin  Y  cos  j 

sin  y  cos  Y 

d 

2 

d 

2 

C2 

1 

kd  Id 

cos  —  cos  — 

2  2 

sinf 

d 

2 

sinf 

d 

2 

D2 

1 

kd  Id 

sin  kd  cos  y 

sin  Id  cos  ^ 

2  2 

d 

d 

A4 

1 

1 

8sin  kd  —  sin  2kd 

6d 

8sin  Id  —  sin  21  d 

6d 

B1 

1 

1 

(-sin2|^  +  27sinf)cosf 

(—  sin  ^  +  27 sin  y )  cos  y 

12d 

\2d 

C4 

1 

kd  Id 

-sin^  +  27sinf 

—  sin  ^  +  27  sin  y 

2  2 

12d 

12d 

1)1 

1 

kd  Id 

(8  sin  kd  —  sin  2kd)  cos  y 

(8  sin  Id  —  sin  2ld)  cos  ^ 

2  2 

6d 

6d 
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Scheme 


Table  III.  Filter  weights  for  finite  elements 
a  =  13  7i 


72 


FET 

FER 


3  +  cos  fcc?  +  2  cos  Y  cos  W  2(sin  A:c/ +  sin  y)  cos /o?  cosysinW 

6  3^  d 

(2  +  cos  kd)(2  +  cos  Id)  (2  +  cos  Ic?)(sin  kd)  (2  +  cos  kd){sm  Id) 
9  Zd  3rf 


The  two  schemes  for  finite  elements  in  Table  III  are  for  triangles  (Finite 
Element-Triangle  or  FET)  and  rectangles  (FER). 

The  Fourier  transform  of  the  steady  state  solution  is  similar  to  the  steady  state 
solution  of  the  continuous  case  and  is  written  as  follows  with  obvious  changes  from 
the  continuous  case. 

Us 

ds 
hs 


\ 

_ 1 

1  ^ 

=  2  2^^ 

A 

Uq 

- 1 

where 


Cd 


gh'yj  -gH'^1'12  -*5'/^72 
-5fi/7i72  gh'yl  igfPli 


(5.8) 


and 

oil'  =  /?/^l  +  Xoili  +  72) 

There  is  a  great  deal  of  similarity  between  the  continuous  case  and  the  discrete  case 
up  to  this  point.  We  can  carry  it  further  and  define  as  before 


D{k,  /,  0)  =  i7i'Uo  +  ^72^0 

C(A;,/,0)  =^72^0- *71^0 

ds;{k,l,0)  =ij'yiho-0vo 

dy{k,l,0)  =  ij'y2ho  +  l3uo- 
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Therefore  the  steady  state  Fourier  transforms  of  the  semi-discrete  shallow  water  equa¬ 
tions  are 


u. 


f2 

=  uo  +  -  (3dy) 


p 

Vs=Vo  +  ^-r(A^i72-D  -  pd^) 


1 

v2  T /2 


(5.9) 


H  p 

hs  =  ho  +  -j-T^{i-1id^  +  i'y^dy). 

Notice  that  the  Rossby  radius  of  deformation  (A)  was  used  instead  of  the  its  discrete 
analog  (A©). 

a.  Finding  the  Steady  State  Solutions  for  Scheme  A2 
Let 

^i{kx+ly) 


For  scheme  A2,  F~^(i‘yiq{k,l))  is 


dkdl. 


=  -f 

2dJ- 


1  p  dq{xo,y) 


d  dx 


dxo 


-  q{-d,y)] 


and  T  ^{p2q{kj))  is 


oo  rd  r 


-  PS 


dqixo,yo)  1  ..  ... 

— d(x  -  xojdyodxo 


ooJ—d  dy  2,d 


-  hi 


1  p  dq(x,yo) 


d  dy 


dyo 


2d 


{q{x,d)-q{x,-d)]. 
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It  follows  in  the  semi-discrete  case  of  the  shallow  water  equations,  by  using  the  con¬ 
volution  theorem  the  steady  state  solution  for  scheme  A2  is 


Usix.y)  =  Uo{x,y)  +  A 


oo  roo  (  -  a,y  -  t)  1 


f  —ooJ  —oo 


dx  2d 


[uo{d,  r,  0)  -  uo(-d,  r,  0)] 


/  1  \  2  I 
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v,(x,y)  =  voix,y)  -|-  A 
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h,{x.y)  =  ho{x,y)  +  ——  /  ^ - - - --[ho(d,r,0)  - /io(-d,T,0)] 
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■l,n{x  -  ff,y-  T)[uo(d,  r,0)  -  uo(-d,r,0)] 


-l,„{x  -  (T,y-  r)[uo(cr,  d,0)  -  uo(a, -d,0)] 
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c)I„,(x  - 


where 


— ^7R*o(c»’,d,0)  —  ho((7,  -d,0)]l  d<Tdr 

oy  f  J 


d  2 

(lx(x,y,0)  =  j-^lho(d,y)  -  hol-d,y)]-Vo 


(5.10) 


Note  that  dx{x,y^Qi)  is  a  centered  difference  approximation  which  approaches  as 
d  approaches  0.  One  should  be  able  to  obtain  similar  results  for  the  other  schemes 
listed  ill  the  tables. 

b.  Solving  for  Im 

The  reason  for  this  thesis  is  to  show  which  semi-discrete  method  has 
an  inverse  Fourier  transform  which  can  be  back  transformed  in  a  closed  form.  By  the 
Lax  equivalence  theorem,  if  we  can  show  that  is  bounded  in  some  sense,  then  the 
solution  is  stable  also.  Numerical  analysis  may  be  used  to  solve  this  problem,  but  a 
simple  closed  formula  is  considered  a  more  ideal  solution.  Let  us  look  again  at 
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dkdl.  (5.11) 


ooJ-ool3^k,  1)  +  X[jUk,  1)  +  7l(^,  01 

Note  that  the  limit  of  Im  for  any  of  the  schemes  when  d  — >■  0  is 

^i{kx+ly) 


/OO  roo 
-OO*/  —OO  1 


-dkdL 


-oo*/-ool  +  A2  (k^  +  P) 

By  using  the  identity  ^  =  /q°°  e~^°'d\,  we  have  the  following  triple  integral, 

Jo  J —ooJ  -‘OO 

which  can  be  rewritten  as 

(  r  H  dg. 

Jo  \*/— OO  J-‘00  J 

Evaluate  the  three  integrals  separately  in  this  form  and  it  easily  can  be  shown  that 

HmA^/,n  =  TpA'o 
d-^0  2tt 

where  Ao  is  the  Bessel  function  of  order  zero  from  (5.6).  This  shows  that  the  steady 
stat#*  solution  of  the  semi-discrete  system  (5.10)  approaches  the  corresponding  solution 
(5.5)  of  the  continuous  case. 

Refering  to  Neta  and  Devito  (see  [Ref,  7])  we  have  sevtTal  second  and 
fourth  order  finite  difference  equations  (Table  II).  In  scheme  A2  we  have 
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+  A2  p 


dkdl. 


(5.12) 
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By  noting  in  (5.12)  that  we  have,  for  example,  cosk  in  the  numerator  and  that  the 
denominator  is  positive  and  finite,  we  see  that  Ia2  does  not  converge.  What  one  can 
do  is  look  at  1^2  to  determine  if  there  is  any  useful  information  which  can  be  drawn 
from  it.  This  information  then  can  be  used  to  predict  some  aspect  of  the  semi-discrete 
shallow  water  equations.  Latta  ([Ref.  31])  suggested  the  following  reasoning.  One 


can  rewrite  Ia2  as 
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G  j—ooj—oo  -ir 
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-f  ^sin^(fcd)  -|-  sin^(W)j 


where  C  =  -rr.  Notice  that  (7  is  a  constant.  Now  if  one  lets  kd  k  and  Id  —y  /,  Ia2 
d^ 

becomes 
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(5.13) 


Let  A;  — >•  A:  -b  TT  and  I  l  +  ir.  This  shows  that  Ia2  is  infinite  where  Ia2  =  e^-^-lA2  and 
Ia2  =  e'^/42.  Therefore  Ia2  is  infinite  when  x  =  2md  and  y  =  2nd  where  m  and  n 
are  integers  and  Ia2  —  0  otherwise.  The  solution  to  Ia2  will  look  like 


Cd^lA2  =  E  E  ^rnnS{k  -  2m)5{l  _  2n)  =  .F  U 


Now  take  the  inverse  Fourier  Transform  of  the  middle  term  to  get 
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V27r/  ^ -b  sin2(A;) -b  sin2(/) 

By  noting  that  we  are  only  interested  in  the  real  part  of  Ia2-,  we  can  write  the  equation 


E  E  (2mA:)  cos  (2n/)  =  ^ 
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Using  the  Fourier  Cosine  series,  we  can  write  the  coefficients  a^n  as 
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Now  using  the  trigonometric  substitution  sin^  /3  = 
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,  we  get 
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which  can  be  written  as 
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Using  the  geometric  series  expansion  ( -  =  l+  x  +  a:^4-a:^  +  ...),  we  can  write 

\  —  X 


OLn 


p2tt  r2'7r  QO  I  C  1  ^ 
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Note  that  cos  (2mA;)  cos  {2nl)  <  1  and  (cos  (2A;)  +  cos  (2/))“"  <  2^.  Therefore 
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This  means  that  a,nn  is  bounded  and  it  follows  that  although  Ia2  infinite,  it  is  stable 
at  each  lattice  point. 


Since  (a  +  b)‘'  =  ^ 
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Using  Gradshteyn  and  Ryshik  (see  [Ref.  32],  p.  374),  one  gets  the  following  definite 
integral. 
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[r  <  2m]; 

[2m  <  r  and  r  —  2m  =  2i]; 

[2m  <  r  and  r  —  2m  =  2i  +  1] 

(5.15) 
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where 

/ 

0  [2m  —  r  =  2i]; 

•s  =  <  1  [2m  —  r  =  +  1]; 

—  1  [2m  —  r  =  4e  —  1] 

Since  amn  is  the  same  across  any  27r  interval,  we  can  write  it  as 


Qmn  =  Xl  (  o  I  o?» )  X]  /  COS  {2km)cos'’ {2k)dk  cos  {2ln)cos''  {2l)dl. 

+  r=o\  r  i 


Noting  that  cosine  is  even,  this  can  be  rewritten  as 


ccmn  =  45^  (2  +  2(7/  ^1  I  Jo  (2A:m)cos’’  (2A:)t;A; cos  (2/n)cos‘'"’'  (20c?/. 

By  using  (5.15),  we  note  that  if  z/  is  odd,  then  amn  =  0.  Also  if  u  is  even  and  r  is  odd 
then  amn  =  0.  Therefore  v  =  (2m  +  2n  +  2i  +  2j)  where  r  =  2m  +  2i.  Now  we  note 
if  r  <  2m  or  1/  —  r  <  2n,  then  amn  =  0.  From  these  observations  of  (5.15),  it  follows 
that 
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This  can  be  reduced  to 
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Where  the  prime  on  the  summation  means  to  count  only  even  indices.  Therefore 
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VI.  CONCLUSIONS 


The  basis  of  this  thesis  was  to  help  the  reader  form  a  better  understanding 
of  the  shallow  water  equations  and  also  to  see  if  an  inverse  Fourier  transform  for 
the  semi-discretized  system  can  give  any  information  on  how  the  system  behaves  as 
a  function  of  time.  The  first  four  chapters  consisted  of  background  material  that  I 
compiled  as  a  result  of  studying  the  shallow  water  equations,  and  also  in  an  attempt  to 
give  the  reader  a  logical  order  to  better  facilitate  understanding  of  the  later  material. 
As  can  be  seen,  various  forms  of  the  shallow  water  equations  were  discussed  and 
derived.  Most  of  that  work  was  taken  directly  from  the  references  and  were  so  noted. 
I  feel  this  thesis  will  stand  as  a  means  for  future  students  to  more  easily  come  to  an 
understanding  of  shallow  water  equations  and  help  them  to  pursue  the  next  step  in 
the  process  by  finding  a  way  of  looking  at  other  finite  difference  and  finite  element 
schemes. 

I  hope  you  enjoyed  reading  this  thesis. 
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